///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/fqmc/rgp_detail.cpp
///
/// \brief   implementation file for some functions used with the restricted geminal power
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<complex>
#include<algorithm>

#include<formic/fqmc/rgp_detail.h>
#include<formic/fqmc/rgp_fortran.h>
#include<formic/numeric/numeric.h>
#include<formic/lapack/interface.h>
#include<formic/exception.h>
#include<formic/ham/ham.h>
#include<formic/timing/timing.h>

#include<formic/output/output.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief computes the 1x1 determinant involved in the single-excitation ratio for the
///        restricted geminal power wavefunction
///
/// \param[in]      n        the number of occupied orbitals
/// \param[in]      d1       the 1st constant, usually equal to zero or one
/// \param[in]      v1       the first v vector, size n
/// \param[in]      u1       the first u vector, size n
/// \param[in]      i1       the first inverse matrix, size n*n
/// \param[out]     w        workspace, size n
///
/// \return  ( d1 + v1^T i1 u1 )
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> S formic::fqmc::rgp_single_ratio(const int n,
                                                    const S d1,
                                                    const S * const v1,
                                                    const S * const u1,
                                                    const S * const i1,
                                                          S * const w)
{

  // i1 u1
  formic::xgemm('T', 'N', n, 1, n, formic::unity(S()), i1, n, u1, n, formic::zero(S()), w, n);

  // d1 + v1^T i1 u1
  const S v1i1u1 = d1 + formic::xdot(n, v1, 1, w, 1);

  return v1i1u1;

}

template double formic::fqmc::rgp_single_ratio(const int n,
                                               const double d1,
                                               const double * const v1,
                                               const double * const u1,
                                               const double * const i1,
                                                     double * const w);
template std::complex<double> formic::fqmc::rgp_single_ratio(const int n,
                                                             const std::complex<double> d1,
                                                             const std::complex<double> * const v1,
                                                             const std::complex<double> * const u1,
                                                             const std::complex<double> * const i1,
                                                                   std::complex<double> * const w);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief computes the 2x2 determinant involved in the double-excitation ratio for the
///        restricted geminal power wavefunction
///
/// \param[in]      n        the number of occupied orbitals
/// \param[in]      d1       the 1st constant, usually equal to zero or one
/// \param[in]      d2       the 2nd constant, usually equal to zero or one
/// \param[in]      v1       the first v vector, size n
/// \param[in]      v2       the second v vector, size n
/// \param[in]      u1       the first u vector, size n
/// \param[in]      u2       the second u vector, size n
/// \param[in]      i1       the first inverse matrix, size n*n
/// \param[in]      i2       the second inverse matrix, size n*n
/// \param[out]     w        workspace, size n
///
/// \return  ( d1 + v1^T i1 u1 ) ( d2 + v2^T i2 u2 ) - ( v2^T i1 u1 ) ( v1^T i2 u2 )
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> S formic::fqmc::rgp_double_ratio(const int n,
                                                    const S d1,
                                                    const S d2,
                                                    const S * const v1,
                                                    const S * const v2,
                                                    const S * const u1,
                                                    const S * const u2,
                                                    const S * const i1,
                                                    const S * const i2,
                                                          S * const w)
{

  // i1 u1
  formic::xgemm('T', 'N', n, 1, n, formic::unity(S()), i1, n, u1, n, formic::zero(S()), w, n);

  // d1 + v1^T i1 u1
  const S v1i1u1 = d1 + formic::xdot(n, v1, 1, w, 1);

  // v2^T i1 u1
  const S v2i1u1 =      formic::xdot(n, v2, 1, w, 1);

  // i2 u2
  formic::xgemm('T', 'N', n, 1, n, formic::unity(S()), i2, n, u2, n, formic::zero(S()), w, n);

  // d2 + v2^T i2 u2
  const S v2i2u2 = d2 + formic::xdot(n, v2, 1, w, 1);

  // v1^T i2 u2
  const S v1i2u2 =      formic::xdot(n, v1, 1, w, 1);

  return v1i1u1 * v2i2u2 - v2i1u1 * v1i2u2;

}

template double formic::fqmc::rgp_double_ratio(const int n,
                                               const double d1,
                                               const double d2,
                                               const double * const v1,
                                               const double * const v2,
                                               const double * const u1,
                                               const double * const u2,
                                               const double * const i1,
                                               const double * const i2,
                                                     double * const w);
template std::complex<double> formic::fqmc::rgp_double_ratio(const int n,
                                                             const std::complex<double> d1,
                                                             const std::complex<double> d2,
                                                             const std::complex<double> * const v1,
                                                             const std::complex<double> * const v2,
                                                             const std::complex<double> * const u1,
                                                             const std::complex<double> * const u2,
                                                             const std::complex<double> * const i1,
                                                             const std::complex<double> * const i2,
                                                                   std::complex<double> * const w);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   builds the occupied pairing matrix
///
/// \param[in]     nap       number of alpha particles
/// \param[in]     nas       number of alpha lattice sites
/// \param[in]     gamma     the all-sites pairing matrix, size nas*nas
/// \param[in]     a_sites   the occupied alpha sites, size nap
/// \param[in]     b_sites   the occupied beta  sites, size nap
/// \param[out]    mat       on exit, the occupied pairing matrix, size nap*nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::fqmc::rgp_build_occupied_matrix(const int nap,
                                                                const int nas,
                                                                const S * const gamma,
                                                                const int * const a_sites,
                                                                const int * const b_sites,
                                                                S * const mat)
{
  for (int i = 0; i < nap; i++)
  for (int j = 0; j < nap; j++)
    mat[i*nap+j] = gamma[ a_sites[i] * nas + b_sites[j] ];
}

template void formic::fqmc::rgp_build_occupied_matrix(const int nap,
                                                      const int nas,
                                                      const double * const gamma,
                                                      const int * const a_sites,
                                                      const int * const b_sites,
                                                      double * const mat);
template void formic::fqmc::rgp_build_occupied_matrix(const int nap,
                                                      const int nas,
                                                      const std::complex<double> * const gamma,
                                                      const int * const a_sites,
                                                      const int * const b_sites,
                                                      std::complex<double> * const mat);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   computes the difference between two rows of the pairing matrix
///
/// \param[in]     nap        number of alpha particles
/// \param[in]     nas        number of alpha lattice sites
/// \param[in]     new_row    the index of the new row, 0 <= new_row < nas
/// \param[in]     old_row    the index of the old row, 0 <= old_row < nas
/// \param[in]     gamma      the all-sites pairing matrix, size nas*nas
/// \param[in]     b_sites    the occupied beta sites, size nap
/// \param[out]    diff       on exit, the difference between the rows, size nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::fqmc::rgp_row_diff(const int nap,
                                                   const int nas,
                                                   const int new_row,
                                                   const int old_row,
                                                   const S * const gamma,
                                                   const int * const b_sites,
                                                   S * const diff)
{
  for (int i = 0; i < nap; i++)
    diff[i] = gamma[ new_row * nas + b_sites[i] ] - gamma[ old_row * nas + b_sites[i] ];
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   computes the difference between two columns of the pairing matrix
///
/// \param[in]     nap        number of alpha particles
/// \param[in]     nas        number of alpha lattice sites
/// \param[in]     new_col    the index of the new column, 0 <= new_col < nas
/// \param[in]     old_col    the index of the old column, 0 <= old_col < nas
/// \param[in]     gamma      the all-sites pairing matrix, size nas*nas
/// \param[in]     a_sites    the occupied alpha sites, size nap
/// \param[out]    diff       on exit, the difference between the columns, size nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::fqmc::rgp_col_diff(const int nap,
                                                   const int nas,
                                                   const int new_col,
                                                   const int old_col,
                                                   const S * const gamma,
                                                   const int * const a_sites,
                                                   S * const diff)
{
  for (int i = 0; i < nap; i++)
    diff[i] = gamma[ a_sites[i] * nas + new_col ] - gamma[ a_sites[i] * nas + old_col ];
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   updates the inverse of the occupied pairing matrix after a single row replacement
///          corresponding to a single alpha particle move
///
/// \param[in]     nap            number of alpha particles
/// \param[in]     nas            number of alpha lattice sites
/// \param[in]     new_occ_site   the newly occupied site, 0 <= new_occ_site < nas
/// \param[in]     old_occ_site   the newly empty site,    0 <= old_occ_site < nas
/// \param[in]     gamma          the all-sites pairing matrix, size nas*nas
/// \param[in]     a_site_to_occ  mapping from alpha sites to their ordering, size nas
/// \param[in]     a_occ_to_site  the occupied alpha sites, size nap
/// \param[in]     b_occ_to_site  the occupied beta sites, size nap
/// \param[in,out] occ_mat        size nap*nap
///                               on entry, the current occupied pairing matrix
///                               on exit, the occupied pairing matrix after the row replacement
/// \param[in,out] inv_mat        size nap*nap
///                               on entry, the current inverse of the occupied pairing matrix
///                               on exit, the inverse after the row replacement
/// \param[out]    wv1            work space, size nap
/// \param[out]    wv2            work space, size nap
/// \param[out]    wv3            work space, size nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::fqmc::rgp_row_update(const int nap,
                                                     const int nas,
                                                     const int new_occ_site,
                                                     const int old_occ_site,
                                                     const S * const gamma,
                                                     const int * const a_site_to_occ,
                                                     const int * const a_occ_to_site,
                                                     const int * const b_occ_to_site,
                                                     S * const occ_mat,
                                                     S * const inv_mat,
                                                     S * const wv1,
                                                     S * const wv2,
                                                     S * const wv3)
{

  // get the index of the row to update
  const int row = a_site_to_occ[old_occ_site];
  assert( row >= 0   );
  assert( row <  nap );

  // get the difference between the new and old row
  formic::fqmc::rgp_row_diff(nap, nas, new_occ_site, old_occ_site, gamma, b_occ_to_site, wv1);

  // update the occupied pairing matrix
  for (int i = 0; i < nap; i++)
    occ_mat[ nap * row + i ] += wv1[i];

  // get the necessary column of the inverse matrix
  for (int i = 0; i < nap; i++)
    wv2[i] = inv_mat[ nap * i + row ];

  // compute the determinant ratio
  const S ratio = formic::unity(S()) + formic::xdot(nap, wv1, 1, wv2, 1);

  // compute the action of the inverse matrix on the row difference
  formic::xgemm('N', 'N', nap, 1, nap, formic::unity(S()), inv_mat, nap, wv1, nap, formic::zero(S()), wv3, nap);

  // update the inverse matrix
  assert( ratio != formic::zero(S()) );
  formic::xgemm('N', 'N', nap, nap, 1, -1.0 / ratio, wv3, nap, wv2, 1, formic::unity(S()), inv_mat, nap);

}

template void formic::fqmc::rgp_row_update(const int nap,
                                           const int nas,
                                           const int new_occ_site,
                                           const int old_occ_site,
                                           const double * const gamma,
                                           const int * const a_site_to_occ,
                                           const int * const a_occ_to_site,
                                           const int * const b_occ_to_site,
                                           double * const occ_mat,
                                           double * const inv_mat,
                                           double * const wv1,
                                           double * const wv2,
                                           double * const wv3);
template void formic::fqmc::rgp_row_update(const int nap,
                                           const int nas,
                                           const int new_occ_site,
                                           const int old_occ_site,
                                           const std::complex<double> * const gamma,
                                           const int * const a_site_to_occ,
                                           const int * const a_occ_to_site,
                                           const int * const b_occ_to_site,
                                           std::complex<double> * const occ_mat,
                                           std::complex<double> * const inv_mat,
                                           std::complex<double> * const wv1,
                                           std::complex<double> * const wv2,
                                           std::complex<double> * const wv3);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   updates the inverse of the occupied pairing matrix after a single column replacement
///          corresponding to a single beta particle move
///
/// \param[in]     nap            number of alpha particles
/// \param[in]     nas            number of alpha lattice sites
/// \param[in]     new_occ_site   the newly occupied site, 0 <= new_occ_site < nas
/// \param[in]     old_occ_site   the newly empty site,    0 <= old_occ_site < nas
/// \param[in]     gamma          the all-sites pairing matrix, size nas*nas
/// \param[in]     b_site_to_occ  mapping from beta sites to their ordering, size nas
/// \param[in]     b_occ_to_site  the occupied beta sites, size nap
/// \param[in]     a_occ_to_site  the occupied alpha sites, size nap
/// \param[in,out] occ_mat        size nap*nap
///                               on entry, the current occupied pairing matrix
///                               on exit, the occupied pairing matrix after the column replacement
/// \param[in,out] inv_mat        size nap*nap
///                               on entry, the current inverse of the occupied pairing matrix
///                               on exit, the inverse after the column replacement
/// \param[out]    wv1            work space, size nap
/// \param[out]    wv2            work space, size nap
/// \param[out]    wv3            work space, size nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void formic::fqmc::rgp_col_update(const int nap,
                                                     const int nas,
                                                     const int new_occ_site,
                                                     const int old_occ_site,
                                                     const S * const gamma,
                                                     const int * const b_site_to_occ,
                                                     const int * const b_occ_to_site,
                                                     const int * const a_occ_to_site,
                                                     S * const occ_mat,
                                                     S * const inv_mat,
                                                     S * const wv1,
                                                     S * const wv2,
                                                     S * const wv3)
{

  // get the index of the column to update
  const int col = b_site_to_occ[old_occ_site];
  assert( col >= 0   );
  assert( col <  nap );

  // get the difference between the new and old column
  formic::fqmc::rgp_col_diff(nap, nas, new_occ_site, old_occ_site, gamma, a_occ_to_site, wv1);

  // update the occupied pairing matrix
  for (int i = 0; i < nap; i++)
    occ_mat[ nap * i + col ] += wv1[i];

  // get the necessary row of the inverse matrix
  for (int i = 0; i < nap; i++)
    wv2[i] = inv_mat[ nap * col + i ];

  // compute the determinant ratio
  const S ratio = formic::unity(S()) + formic::xdot(nap, wv1, 1, wv2, 1);

  // compute the action of the inverse matrix on the column difference
  formic::xgemm('T', 'N', nap, 1, nap, formic::unity(S()), inv_mat, nap, wv1, nap, formic::zero(S()), wv3, nap);

  // update the inverse matrix
  assert( ratio != formic::zero(S()) );
  formic::xgemm('N', 'N', nap, nap, 1, -1.0 / ratio, wv2, nap, wv3, 1, formic::unity(S()), inv_mat, nap);

}

template void formic::fqmc::rgp_col_update(const int nap,
                                           const int nas,
                                           const int new_occ_site,
                                           const int old_occ_site,
                                           const double * const gamma,
                                           const int * const b_site_to_occ,
                                           const int * const b_occ_to_site,
                                           const int * const a_occ_to_site,
                                           double * const occ_mat,
                                           double * const inv_mat,
                                           double * const wv1,
                                           double * const wv2,
                                           double * const wv3);
template void formic::fqmc::rgp_col_update(const int nap,
                                           const int nas,
                                           const int new_occ_site,
                                           const int old_occ_site,
                                           const std::complex<double> * const gamma,
                                           const int * const b_site_to_occ,
                                           const int * const b_occ_to_site,
                                           const int * const a_occ_to_site,
                                           std::complex<double> * const occ_mat,
                                           std::complex<double> * const inv_mat,
                                           std::complex<double> * const wv1,
                                           std::complex<double> * const wv2,
                                           std::complex<double> * const wv3);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   returns the ratio between the AGP determinant for a new configuration and the
///          current configuration
///
/// \param[in]     nap            number of alpha particles
/// \param[in]     nas            number of alpha lattice sites
/// \param[in]     nds            the number of sites that differ between the configurations
/// \param[in]     newly_occ      newly occupied lattice sites, size nds/2
/// \param[in]     newly_emp      newly empty lattice sites, size nds/2
/// \param[in]     gamma          the all-sites pairing matrix, size nas*nas
/// \param[in]     inv_mat        inverse of the occupied pairing matrix, size nap*nap
/// \param[in]     der_coeff      derivative coefficients for pairing matrix variables, size nas*nas
/// \param[in]     cont_inv       contraction of the inverse matrix with the der. coeff., size nap*nap
/// \param[in]     a_occ_to_site  the occupied alpha sites, size nap
/// \param[in]     b_occ_to_site  the occupied beta sites, size nap
/// \param[in]     a_site_to_occ  mapping from alpha sites to their ordering, size nas
/// \param[in]     b_site_to_occ  mapping from beta sites to their ordering, size nas
/// \param[out]    work           work space, size 5*nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> S formic::fqmc::rgp_det_ratio_deriv(const int nap,
                                                       const int nas,
                                                       const int nds,
                                                       const int * const newly_occ,
                                                       const int * const newly_emp,
                                                       const S * const gamma,
                                                       const S * const inv_mat,
                                                       const S * const der_coeff,
                                                       const S * const cont_inv,
                                                       const int * const a_occ_to_site,
                                                       const int * const b_occ_to_site,
                                                       const int * const a_site_to_occ,
                                                       const int * const b_site_to_occ,
                                                       S * const work)
{

  // get workspace vectors
  S * const wv1 = work + 0*nap;
  S * const wv2 = work + 1*nap;
  S * const wv3 = work + 2*nap;
  S * const wv4 = work + 3*nap;
  S * const wv5 = work + 4*nap;

  // check that the number of sites affected by the move is supported
  if ( nds != 0 && nds != 2 && nds != 4 )
    throw formic::Exception("unsupported number of differing sites nds = %i in formic::fqmc::rgp_det_ratio_deriv") % nds;

  // check whether we are computing the contraction of the differential ratios with derivative coefficients or just the ratio
  const bool do_der_cont = ( der_coeff != 0 && cont_inv != 0 );
  //if ( do_der_cont )
  //  throw formic::Exception("derivative contraction not yet implemented in formic::fqmc::rgp_det_ratio_deriv");

  // initialize a variable that will hold the derivative contraction term corresponding to the derivative hitting the regular determinant
  S dc_reg_det = formic::zero(S());
  if ( do_der_cont )
    for ( int i = 0; i < nap; i++)
      for ( int j = 0; j < nap; j++)
        dc_reg_det += inv_mat[nap*j+i] * der_coeff[ a_occ_to_site[i] * nas + b_occ_to_site[j] ];

  // special handling when no sites are changed
  if ( nds == 0 && do_der_cont )
    return dc_reg_det;
  else if ( nds == 0 )
    return formic::unity(S());

  // determine which among the affected lattice sites are newly occupied, and which were previously occupied
  int n_a_diff = 0;
  int n_b_diff = 0;
  int a_new_occ_sites[2];
  int b_new_occ_sites[2];
  int a_old_occ_sites[2];
  int b_old_occ_sites[2];
  for (int i = 0; 2*i < nds; i++) {
    if ( newly_occ[i] % 2 != newly_emp[i] % 2 )
      throw formic::Exception("spin mismatch in occ/empty pair number %i in formic::fqmc::rgp_det_ratio_deriv") % i;
    if ( newly_occ[i] % 2 == 0 ) {
      a_new_occ_sites[n_a_diff] = newly_occ[i] / 2;
      a_old_occ_sites[n_a_diff] = newly_emp[i] / 2;
      n_a_diff++;
    } else {
      b_new_occ_sites[n_b_diff] = newly_occ[i] / 2;
      b_old_occ_sites[n_b_diff] = newly_emp[i] / 2;
      n_b_diff++;
    }
  }
  n_a_diff *= 2;
  n_b_diff *= 2;

  // create the return value
  S retval = formic::zero(S());

  // single alpha move
  if ( n_a_diff == 2 && n_b_diff == 0 ) {

    // determine which column of the inverse matrix we need
    const int col = a_site_to_occ[a_old_occ_sites[0]];

    // create the column vector
    std::fill(wv2, wv2 + nap, formic::zero(S()));
    wv2[col] = formic::unity(S());

    // get the difference between the new and the old rows
    formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_occ_to_site, wv1);

    // derivative contraction
    if ( do_der_cont ) {

      //// derivative hits the regular determinant
      //retval += dc_reg_det * formic::fqmc::rgp_single_ratio(nap, formic::unity(S()), wv1, wv2, inv_mat,  wv3);

      // derivative hits the inverse matrix
      retval +=        1.0 * formic::fqmc::rgp_single_ratio(nap, formic::zero(S()),  wv1, wv2, cont_inv, wv3);

      // derivative hits the row difference vector
      formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], der_coeff, b_occ_to_site, wv1);
      retval +=        1.0 * formic::fqmc::rgp_single_ratio(nap, formic::zero(S()),  wv1, wv2, inv_mat,  wv3);

    // just the ratio
    } else {

      // compute the ratio of the new and old determinants
      retval = formic::fqmc::rgp_single_ratio(nap, formic::unity(S()), wv1, wv2, inv_mat, wv3);

      //// compute the ratio of the new and old determinants by multiplying the row difference by the corresponding column of the inverse matrix
      //retval = formic::unity(S()) + formic::xdot(nap, inv_mat + col, nap, wv1, 1);

    }

  // single beta move
  } else if ( n_a_diff == 0 && n_b_diff == 2 ) {

    // determine which row of the inverse matrix we need
    const int row = b_site_to_occ[b_old_occ_sites[0]];

    // create the row vector
    std::fill(wv2, wv2 + nap, formic::zero(S()));
    wv2[row] = formic::unity(S());

    // get the difference between the new and the old columns
    formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_occ_to_site, wv1);

    // derivative contraction
    if ( do_der_cont ) {

      //// derivative hits the regular determinant
      //retval += dc_reg_det * formic::fqmc::rgp_single_ratio(nap, formic::unity(S()), wv2, wv1,  inv_mat, wv3);

      // derivative hits the inverse matrix
      retval +=        1.0 * formic::fqmc::rgp_single_ratio(nap, formic::zero(S()),  wv2, wv1, cont_inv, wv3);

      // derivative hits the column difference vector
      formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], der_coeff, a_occ_to_site, wv1);
      retval +=        1.0 * formic::fqmc::rgp_single_ratio(nap, formic::zero(S()),  wv2, wv1,  inv_mat, wv3);

    // just the ratio
    } else {

      // compute the ratio of the new and old determinants
      retval = formic::fqmc::rgp_single_ratio(nap, formic::unity(S()), wv2, wv1, inv_mat, wv3);

      //// compute the ratio of the new and old determinants by multiplying the column difference by the corresponding row of the inverse matrix
      //retval = formic::unity(S()) + formic::xdot(nap, inv_mat + nap*row, 1, wv1, 1);

    }

  // double alpha move
  } else if ( n_a_diff == 4 && n_b_diff == 0 ) {

    // determine which columns of the inverse matrix we need
    const int col1 = a_site_to_occ[a_old_occ_sites[0]];
    const int col2 = a_site_to_occ[a_old_occ_sites[1]];
    assert( col1 != col2 );

    // create the column vectors
    std::fill(wv3, wv3 + nap, formic::zero(S()));
    std::fill(wv4, wv4 + nap, formic::zero(S()));
    wv3[col1] = formic::unity(S());
    wv4[col2] = formic::unity(S());

    // get the differences between the new and the old rows
    formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_occ_to_site, wv1);
    formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[1], a_old_occ_sites[1], gamma, b_occ_to_site, wv2);

    // derivative contraction
    if ( do_der_cont ) {

      //// derivative hits the regular determinant
      //retval += dc_reg_det * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::unity(S()), wv1, wv2, wv3, wv4, inv_mat,  inv_mat,  wv5);

      // derivative hits the 1st inverse matrix
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::zero(S()),  formic::unity(S()), wv1, wv2, wv3, wv4, cont_inv, inv_mat,  wv5);

      // derivative hits the 2nd inverse matrix
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::zero(S()),  wv1, wv2, wv3, wv4, inv_mat,  cont_inv, wv5);

      // derivative hits the 1st row difference vector
      formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], der_coeff, b_occ_to_site, wv1);
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::zero(S()),  formic::unity(S()), wv1, wv2, wv3, wv4, inv_mat,  inv_mat,  wv5);

      // derivative hits the 2nd row difference vector
      formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma,     b_occ_to_site, wv1);
      formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[1], a_old_occ_sites[1], der_coeff, b_occ_to_site, wv2);
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::zero(S()),  wv1, wv2, wv3, wv4, inv_mat,  inv_mat,  wv5);

    // just the ratio
    } else {

      // compute the ratio of the new and old determinants
      retval = formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::unity(S()), wv1, wv2, wv3, wv4, inv_mat, inv_mat, wv5);
      //retval =    ( formic::unity(S()) + formic::xdot(nap, inv_mat + col1, nap, wv1, 1) )
      //          * ( formic::unity(S()) + formic::xdot(nap, inv_mat + col2, nap, wv2, 1) )
      //          - formic::xdot(nap, inv_mat + col1, nap, wv2, 1) * formic::xdot(nap, inv_mat + col2, nap, wv1, 1);

    }

  // double beta move
  } else if ( n_a_diff == 0 && n_b_diff == 4 ) {

    // determine which row of the inverse matrix we need
    const int row1 = b_site_to_occ[b_old_occ_sites[0]];
    const int row2 = b_site_to_occ[b_old_occ_sites[1]];
    assert( row1 != row2 );

    // create the row vectors
    std::fill(wv3, wv3 + nap, formic::zero(S()));
    std::fill(wv4, wv4 + nap, formic::zero(S()));
    wv3[row1] = formic::unity(S());
    wv4[row2] = formic::unity(S());

    // get the differences between the new and the old columns
    formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_occ_to_site, wv1);
    formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[1], b_old_occ_sites[1], gamma, a_occ_to_site, wv2);

    // derivative contraction
    if ( do_der_cont ) {

      //// derivative hits the regular determinant
      //retval += dc_reg_det * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::unity(S()), wv3, wv4, wv1, wv2, inv_mat,  inv_mat,  wv5);

      // derivative hits the 1st inverse matrix
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::zero(S()),  formic::unity(S()), wv3, wv4, wv1, wv2, cont_inv, inv_mat,  wv5);

      // derivative hits the 2nd inverse matrix
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::zero(S()),  wv3, wv4, wv1, wv2, inv_mat,  cont_inv, wv5);

      // derivative hits the 1st column difference vector
      formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], der_coeff, a_occ_to_site, wv1);
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::zero(S()),  formic::unity(S()), wv3, wv4, wv1, wv2, inv_mat,  inv_mat,  wv5);

      // derivative hits the 2nd column difference vector
      formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma,     a_occ_to_site, wv1);
      formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[1], b_old_occ_sites[1], der_coeff, a_occ_to_site, wv2);
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::zero(S()),  wv3, wv4, wv1, wv2, inv_mat,  inv_mat,  wv5);

    // just the ratio
    } else {

      // compute the ratio of the new and old determinants
      retval = formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::unity(S()), wv3, wv4, wv1, wv2, inv_mat, inv_mat, wv5);
      //retval =    ( formic::unity(S()) + formic::xdot(nap, inv_mat + nap*row1, 1, wv1, 1) )
      //          * ( formic::unity(S()) + formic::xdot(nap, inv_mat + nap*row2, 1, wv2, 1) )
      //          - formic::xdot(nap, inv_mat + nap*row1, 1, wv2, 1) * formic::xdot(nap, inv_mat + nap*row2, 1, wv1, 1);

    }

  // alpha beta move
  } else if ( n_a_diff == 2 && n_b_diff == 2 ) {

    // determine which column and row of the inverse matrix we need
    const int col = a_site_to_occ[a_old_occ_sites[0]];
    const int row = b_site_to_occ[b_old_occ_sites[0]];

    // create the unit row/column vectors
    std::fill(wv3, wv3 + nap, formic::zero(S()));
    std::fill(wv4, wv4 + nap, formic::zero(S()));
    wv3[col] = formic::unity(S());
    wv4[row] = formic::unity(S());

    // get the differences between the new and the old rows
    formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_occ_to_site, wv1);

    // get the differences between the new and the old columns
    formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_occ_to_site, wv2);

    // special handling of the element at which the row and column cross
    wv1[row] = formic::zero(S());
    wv2[col] = gamma[ a_new_occ_sites[0] * nas + b_new_occ_sites[0] ] - gamma[ a_old_occ_sites[0] * nas + b_old_occ_sites[0] ];

    // derivative contraction
    if ( do_der_cont ) {

      //// derivative hits the regular determinant
      //retval += dc_reg_det * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::unity(S()), wv1, wv4, wv3, wv2, inv_mat,  inv_mat,  wv5);

      // derivative hits the 1st inverse matrix
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::zero(S()),  formic::unity(S()), wv1, wv4, wv3, wv2, cont_inv, inv_mat,  wv5);

      // derivative hits the 2nd inverse matrix
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::zero(S()),  wv1, wv4, wv3, wv2, inv_mat,  cont_inv, wv5);

      // derivative hits the row difference vector
      formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], der_coeff, b_occ_to_site, wv1);
      wv1[row] = formic::zero(S());
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::zero(S()),  formic::unity(S()), wv1, wv4, wv3, wv2, inv_mat,  inv_mat,  wv5);

      // derivative hits the column difference vector
      formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma,     b_occ_to_site, wv1);
      formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], der_coeff, a_occ_to_site, wv2);
      wv1[row] = formic::zero(S());
      wv2[col] = der_coeff[ a_new_occ_sites[0] * nas + b_new_occ_sites[0] ] - der_coeff[ a_old_occ_sites[0] * nas + b_old_occ_sites[0] ];
      retval +=        1.0 * formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::zero(S()),  wv1, wv4, wv3, wv2, inv_mat,  inv_mat,  wv5);

    // just the ratio
    } else {

      // compute the ratio of the new and old determinants
      retval = formic::fqmc::rgp_double_ratio(nap, formic::unity(S()), formic::unity(S()), wv1, wv4, wv3, wv2, inv_mat, inv_mat, wv5);

      //// compute the action of the inverse matrix on the column difference
      //formic::xgemm('T', 'N', nap, 1, nap, formic::unity(S()), inv_mat, nap, wv2, nap, formic::zero(S()), wv3, nap);

      //// compute the ratio of the new and old determinants
      //retval =   ( formic::unity(S()) + formic::xdot(nap, inv_mat + col, nap, wv1, 1) )
      //         * ( formic::unity(S()) + formic::xdot(nap, inv_mat + nap*row, 1, wv2, 1) )
      //         - inv_mat[ nap * row + col ] * formic::xdot(nap, wv1, 1, wv3, 1);

    }

//    // check that the alpha and beta operators are arranged in one of the four allowed orders
//    const int * const ds = diff_sites;
//    assert( ( ds[0] % 2 == 0 && ds[1] % 2 == 0 && ds[2] % 2 == 1 && ds[3] % 2 == 1 ) ||    // RLRL = aabb
//            ( ds[0] % 2 == 1 && ds[1] % 2 == 1 && ds[2] % 2 == 0 && ds[3] % 2 == 0 ) ||    // RLRL = bbaa
//            ( ds[0] % 2 == 0 && ds[1] % 2 == 1 && ds[2] % 2 == 1 && ds[3] % 2 == 0 ) ||    // RLRL = abba
//            ( ds[0] % 2 == 1 && ds[1] % 2 == 0 && ds[2] % 2 == 0 && ds[3] % 2 == 1 )    ); // RLRL = baab
//
//    // if the operators are in the abba or baab ordering, apply a sign since we have assumed aabb ordering
//    if ( ds[0] % 2 != ds[1] % 2 )
//      retval = -retval;

  // unknown move type
  } else throw formic::Exception("unknown move type in formic::fqmc::rgp_det_ratio_deriv with nds = %i") % nds;

  // return the determinant ratio
  return retval;

}

template double formic::fqmc::rgp_det_ratio_deriv(const int nap,
                                                  const int nas,
                                                  const int nds,
                                                  const int * const newly_occ,
                                                  const int * const newly_emp,
                                                  const double * const gamma,
                                                  const double * const inv_mat,
                                                  const double * const der_coeff,
                                                  const double * const cont_inv,
                                                  const int * const a_occ_to_site,
                                                  const int * const b_occ_to_site,
                                                  const int * const a_site_to_occ,
                                                  const int * const b_site_to_occ,
                                                  double * const work);
template std::complex<double> formic::fqmc::rgp_det_ratio_deriv(const int nap,
                                                                const int nas,
                                                                const int nds,
                                                                const int * const newly_occ,
                                                                const int * const newly_emp,
                                                                const std::complex<double> * const gamma,
                                                                const std::complex<double> * const inv_mat,
                                                                const std::complex<double> * const der_coeff,
                                                                const std::complex<double> * const cont_inv,
                                                                const int * const a_occ_to_site,
                                                                const int * const b_occ_to_site,
                                                                const int * const a_site_to_occ,
                                                                const int * const b_site_to_occ,
                                                                std::complex<double> * const work);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   returns the ratio between the AGP determinant for a new configuration and the
///          current configuration
///
/// \param[in]     nap            number of alpha particles
/// \param[in]     nas            number of alpha lattice sites
/// \param[in]     nds            the number of sites that differ between the configurations
/// \param[in]     newly_occ      newly occupied lattice sites, size nds/2
/// \param[in]     newly_emp      newly empty lattice sites, size nds/2
/// \param[in]     gamma          the all-sites pairing matrix, size nas*nas
/// \param[in]     inv_mat        inverse of the occupied pairing matrix, size nap*nap
/// \param[in]     a_occ_to_site  the occupied alpha sites, size nap
/// \param[in]     b_occ_to_site  the occupied beta sites, size nap
/// \param[in]     a_site_to_occ  mapping from alpha sites to their ordering, size nas
/// \param[in]     b_site_to_occ  mapping from beta sites to their ordering, size nas
/// \param[out]    wv1            work space, size nap
/// \param[out]    wv2            work space, size nap
/// \param[out]    wv3            work space, size nap
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> S formic::fqmc::rgp_det_ratio(const int nap,
                                                 const int nas,
                                                 const int nds,
                                                 const int * const newly_occ,
                                                 const int * const newly_emp,
                                                 const S * const gamma,
                                                 const S * const inv_mat,
                                                 const int * const a_occ_to_site,
                                                 const int * const b_occ_to_site,
                                                 const int * const a_site_to_occ,
                                                 const int * const b_site_to_occ,
                                                 S * const wv1,
                                                 S * const wv2,
                                                 S * const wv3)
{

  // check that the number of sites affected by the move is supported
  if ( nds != 0 && nds != 2 && nds != 4 )
    throw formic::Exception("unsupported number of differing sites nds = %i in formic::fqmc::rgp_det_ratio") % nds;

  // return one if no sites were changed
  if (nds == 0)
    return formic::unity(S());

  // determine which among the affected lattice sites are newly occupied, and which were previously occupied
  int n_a_diff = 0;
  int n_b_diff = 0;
  int a_new_occ_sites[2];
  int b_new_occ_sites[2];
  int a_old_occ_sites[2];
  int b_old_occ_sites[2];
  for (int i = 0; 2*i < nds; i++) {
    if ( newly_occ[i] % 2 != newly_emp[i] % 2 )
      throw formic::Exception("spin mismatch in occ/empty pair number %i in formic::fqmc::rgp_det_ratio") % i;
    if ( newly_occ[i] % 2 == 0 ) {
      a_new_occ_sites[n_a_diff] = newly_occ[i] / 2;
      a_old_occ_sites[n_a_diff] = newly_emp[i] / 2;
      n_a_diff++;
    } else {
      b_new_occ_sites[n_b_diff] = newly_occ[i] / 2;
      b_old_occ_sites[n_b_diff] = newly_emp[i] / 2;
      n_b_diff++;
    }
  }
  n_a_diff *= 2;
  n_b_diff *= 2;

  // create the return value
  S retval;

  // single alpha move
  if ( n_a_diff == 2 && n_b_diff == 0 ) {

    // get the difference between the new and the old rows
    formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_occ_to_site, wv1);

    // determine which column of the inverse matrix we need
    const int col = a_site_to_occ[a_old_occ_sites[0]];

    // compute the ratio of the new and old determinants by multiplying the row difference by the corresponding column of the inverse matrix
    retval = formic::unity(S()) + formic::xdot(nap, inv_mat + col, nap, wv1, 1);

  // single beta move
  } else if ( n_a_diff == 0 && n_b_diff == 2 ) {

    // get the difference between the new and the old columns
    formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_occ_to_site, wv1);

    // determine which row of the inverse matrix we need
    const int row = b_site_to_occ[b_old_occ_sites[0]];

    // compute the ratio of the new and old determinants by multiplying the column difference by the corresponding row of the inverse matrix
    retval = formic::unity(S()) + formic::xdot(nap, inv_mat + nap*row, 1, wv1, 1);

  // double alpha move
  } else if ( n_a_diff == 4 && n_b_diff == 0 ) {

    // get the differences between the new and the old rows
    formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_occ_to_site, wv1);
    formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[1], a_old_occ_sites[1], gamma, b_occ_to_site, wv2);

    // determine which columns of the inverse matrix we need
    const int col1 = a_site_to_occ[a_old_occ_sites[0]];
    const int col2 = a_site_to_occ[a_old_occ_sites[1]];
    assert( col1 != col2 );

    // compute the ratio of the new and old determinants
    retval =    ( formic::unity(S()) + formic::xdot(nap, inv_mat + col1, nap, wv1, 1) )
              * ( formic::unity(S()) + formic::xdot(nap, inv_mat + col2, nap, wv2, 1) )
              - formic::xdot(nap, inv_mat + col1, nap, wv2, 1) * formic::xdot(nap, inv_mat + col2, nap, wv1, 1);

  // double beta move
  } else if ( n_a_diff == 0 && n_b_diff == 4 ) {

    // get the differences between the new and the old columns
    formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_occ_to_site, wv1);
    formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[1], b_old_occ_sites[1], gamma, a_occ_to_site, wv2);

    // determine which row of the inverse matrix we need
    const int row1 = b_site_to_occ[b_old_occ_sites[0]];
    const int row2 = b_site_to_occ[b_old_occ_sites[1]];
    assert( row1 != row2 );

    // compute the ratio of the new and old determinants
    retval =    ( formic::unity(S()) + formic::xdot(nap, inv_mat + nap*row1, 1, wv1, 1) )
              * ( formic::unity(S()) + formic::xdot(nap, inv_mat + nap*row2, 1, wv2, 1) )
              - formic::xdot(nap, inv_mat + nap*row1, 1, wv2, 1) * formic::xdot(nap, inv_mat + nap*row2, 1, wv1, 1);

  // alpha beta move
  } else if ( n_a_diff == 2 && n_b_diff == 2 ) {

    // get the differences between the new and the old rows
    formic::fqmc::rgp_row_diff(nap, nas, a_new_occ_sites[0], a_old_occ_sites[0], gamma, b_occ_to_site, wv1);

    // get the differences between the new and the old columns
    formic::fqmc::rgp_col_diff(nap, nas, b_new_occ_sites[0], b_old_occ_sites[0], gamma, a_occ_to_site, wv2);

    // determine which column and row of the inverse matrix we need
    const int col = a_site_to_occ[a_old_occ_sites[0]];
    const int row = b_site_to_occ[b_old_occ_sites[0]];

    // special handling of the element at which the row and column cross
    wv1[row] = formic::zero(S());
    wv2[col] = gamma[ a_new_occ_sites[0] * nas + b_new_occ_sites[0] ] - gamma[ a_old_occ_sites[0] * nas + b_old_occ_sites[0] ];

    // compute the action of the inverse matrix on the column difference
    formic::xgemm('T', 'N', nap, 1, nap, formic::unity(S()), inv_mat, nap, wv2, nap, formic::zero(S()), wv3, nap);

    // compute the ratio of the new and old determinants
    retval =   ( formic::unity(S()) + formic::xdot(nap, inv_mat + col, nap, wv1, 1) )
             * ( formic::unity(S()) + formic::xdot(nap, inv_mat + nap*row, 1, wv2, 1) )
             - inv_mat[ nap * row + col ] * formic::xdot(nap, wv1, 1, wv3, 1);

//    // check that the alpha and beta operators are arranged in one of the four allowed orders
//    const int * const ds = diff_sites;
//    assert( ( ds[0] % 2 == 0 && ds[1] % 2 == 0 && ds[2] % 2 == 1 && ds[3] % 2 == 1 ) ||    // RLRL = aabb
//            ( ds[0] % 2 == 1 && ds[1] % 2 == 1 && ds[2] % 2 == 0 && ds[3] % 2 == 0 ) ||    // RLRL = bbaa
//            ( ds[0] % 2 == 0 && ds[1] % 2 == 1 && ds[2] % 2 == 1 && ds[3] % 2 == 0 ) ||    // RLRL = abba
//            ( ds[0] % 2 == 1 && ds[1] % 2 == 0 && ds[2] % 2 == 0 && ds[3] % 2 == 1 )    ); // RLRL = baab
//
//    // if the operators are in the abba or baab ordering, apply a sign since we have assumed aabb ordering
//    if ( ds[0] % 2 != ds[1] % 2 )
//      retval = -retval;

  // unknown move type
  } else throw formic::Exception("unknown move type in formic::fqmc::rgp_det_ratio with nds = %i") % nds;

  // return the determinant ratio
  return retval;

}

template double formic::fqmc::rgp_det_ratio(const int nap,
                                       const int nas,
                                       const int nds,
                                       const int * const newly_occ,
                                       const int * const newly_emp,
                                       const double * const gamma,
                                       const double * const inv_mat,
                                       const int * const a_occ_to_site,
                                       const int * const b_occ_to_site,
                                       const int * const a_site_to_occ,
                                       const int * const b_site_to_occ,
                                       double * const wv1,
                                       double * const wv2,
                                       double * const wv3);
template std::complex<double> formic::fqmc::rgp_det_ratio(const int nap,
                                                          const int nas,
                                                          const int nds,
                                                          const int * const newly_occ,
                                                          const int * const newly_emp,
                                                          const std::complex<double> * const gamma,
                                                          const std::complex<double> * const inv_mat,
                                                          const int * const a_occ_to_site,
                                                          const int * const b_occ_to_site,
                                                          const int * const a_site_to_occ,
                                                          const int * const b_site_to_occ,
                                                          std::complex<double> * const wv1,
                                                          std::complex<double> * const wv2,
                                                          std::complex<double> * const wv3);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   returns the local energy (or the local hamiltonian derivative contraction)
///          for the initialized configuration
///
/// \param[in]     nap              number of alpha particles
/// \param[in]     nas              number of alpha lattice sites
/// \param[in]     ham              the Hamiltonian
/// \param[in]     gamma            the all-sites pairing matrix, size nas*nas
/// \param[in]     inv_mat          inverse of the occupied pairing matrix, size nap*nap
/// \param[in]     der_coeff        derivative coefficients for pairing matrix variables, size nas*nas
/// \param[in]     cont_inv         contraction of the inverse matrix with the der. coeff., size nap*nap
/// \param[in]     a_occ_to_site    the occupied alpha sites, size nap
/// \param[in]     b_occ_to_site    the occupied beta sites, size nap
/// \param[in]     a_site_to_occ    mapping from alpha sites to their ordering, size nas
/// \param[in]     b_site_to_occ    mapping from beta sites to their ordering, size nas
/// \param[in]     a_unocc_to_site  the unoccupied alpha sites, size nas-nap
/// \param[in]     b_unocc_to_site  the unoccupied beta sites, size nas-nap
/// \param[in]     a_site_to_unocc  mapping from alpha sites to their ordering, size nas
/// \param[in]     b_site_to_unocc  mapping from beta sites to their ordering, size nas
/// \param[out]    work             work space, size 5*nap
///
///  If either der_coeff or cont_inv are null, then the local energy is computed.
///  Otherwise, the Hamiltonian-based derivative ratio contraction using these derivatives
///  and their contracted inverse matrix is computed.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> S formic::fqmc::rgp_simple_local_energy(const int nap,
                                                           const int nas,
                                                           const formic::Hamiltonian<S> & ham,
                                                           const S * const gamma,
                                                           const S * const inv_mat,
                                                           const S * const der_coeff,
                                                           const S * const cont_inv,
                                                           const int * const a_occ_to_site,
                                                           const int * const b_occ_to_site,
                                                           const int * const a_site_to_occ,
                                                           const int * const b_site_to_occ,
                                                           const int * const a_unocc_to_site,
                                                           const int * const b_unocc_to_site,
                                                           const int * const a_site_to_unocc,
                                                           const int * const b_site_to_unocc,
                                                           S * const work)
{

  // local struct to hold some local functions
  struct local_funcs {

//    static void prep_aaaa(
//                          int no,
//                          int nv,
//                          double * tia,
//                          double * ti,
//                          double * viajb,
//                          double * vijb,
//                          double * vij,
//                          double * via,
//                          double * vi,
//                          double * Vmat,
//                          double * Wmat,
//                          double * Bmat,
//                          double * VBmat,
//                          double * WBmat)
//    {
//      FORMIC_FC_GLOBAL_(rgp_fort_energy_prep_aaaa,RGP_FORT_ENERGY_PREP_AAAA)
//        (&no, &nv, tia, ti, viajb, vijb, vij, via, vi, Vmat, Wmat, Bmat, VBmat, WBmat);
//    }

//    static void prep_aaaa(
//                          int no,
//                          int nv,
//                          std::complex<double> * tia,
//                          std::complex<double> * ti,
//                          std::complex<double> * viajb,
//                          std::complex<double> * vijb,
//                          std::complex<double> * vij,
//                          std::complex<double> * via,
//                          std::complex<double> * vi,
//                          std::complex<double> * Vmat,
//                          std::complex<double> * Wmat,
//                          std::complex<double> * Bmat,
//                          std::complex<double> * VBmat,
//                          std::complex<double> * WBmat)
//    {
//      throw formic::Exception("no implementation of fast rgp hamiltonian contractions for complex numbers");
//    }

    static void get_aaaa(
                         int no,
                         int nv,
                         double * tia,
                         //double * ti,
                         double * viajb,
                         //double * vijb,
                         //double * vij,
                         //double * via,
                         //double * vi,
                         double * Bmat,
                         double * Vmat,
                         double * VBmat,
                         //double * WBmat,
                         //double * wia,
                         //double * wij,
                         double * waj,
                         double * retval)
    {
      throw formic::Exception("no longer using get_aaaa for energy");
      //FORMIC_FC_GLOBAL_(rgp_fortran_energy_aaaa,RGP_FORTRAN_ENERGY_AAAA)
      //  (&no, &nv, tia, viajb, Bmat, Vmat, VBmat, waj, retval);
      //  //(&no, &nv, tia, ti, viajb, vijb, vij, via, vi, VBmat, WBmat, wia, wij, waj, retval);
    }

    static void get_aaaa(
                         int no,
                         int nv,
                         std::complex<double> * tia,
                         //std::complex<double> * ti,
                         std::complex<double> * viajb,
                         //std::complex<double> * vijb,
                         //std::complex<double> * vij,
                         //std::complex<double> * via,
                         //std::complex<double> * vi,
                         std::complex<double> * Bmat,
                         std::complex<double> * Vmat,
                         std::complex<double> * VBmat,
                         //std::complex<double> * WBmat,
                         //std::complex<double> * wia,
                         //std::complex<double> * wij,
                         std::complex<double> * waj,
                         std::complex<double> * retval)
    {
      throw formic::Exception("no implementation of fast rgp hamiltonian contractions for complex numbers");
    }

    static void get_aabb(
                         int no,
                         int nv,
                 double * Tiajb,
                double * newrow,
                double * newcol,
                  double * Xmat,
                  double * Bmat,
                  double * Btrn,
                  double * VBia,
                  double * BUkb,
                 double * VBUab,
                  double * TBab,
                  double * retval)
    {
      throw formic::Exception("no longer using get_aabb for energy");
      //FORMIC_FC_GLOBAL_(rgp_fortran_energy_aabb,RGP_FORTRAN_ENERGY_AABB)
      //  (&no, &nv, Tiajb, newrow, newcol, Xmat, Bmat, Btrn, VBia, BUkb, VBUab, TBab, retval);
    }

    static void get_aabb(
                         int no,
                         int nv,
                 std::complex<double> * Tiajb,
                std::complex<double> * newrow,
                std::complex<double> * newcol,
                  std::complex<double> * Xmat,
                  std::complex<double> * Bmat,
                  std::complex<double> * Btrn,
                  std::complex<double> * VBia,
                  std::complex<double> * BUkb,
                 std::complex<double> * VBUab,
                  std::complex<double> * TBab,
                  std::complex<double> * retval)
    {
      throw formic::Exception("no implementation of fast rgp hamiltonian contractions for complex numbers");
    }

    // function to compute energy from same-spin Hamiltonian operators
    static S same_spin(const int beta,
                       const int nap,
                       const int nas,
                       const formic::Hamiltonian<S> & ham,
                       const S * const gamma,
                       const S * const inv_mat,
                       const S * const der_coeff,
                       const S * const cont_inv,
                       const int * const occ_to_site,
                       const int * const unocc_to_site,
                       const int * const a_occ_to_site,
                       const int * const b_occ_to_site,
                       const int * const a_site_to_occ,
                       const int * const b_site_to_occ,
                       S * const work)
    {

      const int nau = nas - nap;

      // get work vectors
      S * const wv1 = work + 0*nap;
      S * const wv2 = work + 1*nap;
      S * const wv3 = work + 2*nap;

      // initialize contribution
      S retval = formic::zero(S());

      // get integer arrays describe the excitations
      int newly_full[] = {0,0};
      int newly_empt[] = {0,0};

      // diagonal one and two electron terms
      if ( der_coeff == 0 || cont_inv == 0 ) {

        const S rat = formic::fqmc::rgp_det_ratio_deriv(nap,
                                                        nas,
                                                        0,
                                                        newly_full,
                                                        newly_empt,
                                                        gamma,
                                                        inv_mat,
                                                        der_coeff,
                                                        cont_inv,
                                                        a_occ_to_site,
                                                        b_occ_to_site,
                                                        a_site_to_occ,
                                                        b_site_to_occ,
                                                        work);

        for (int i = 0; i < nap; i++) {

          // ii
          retval += rat * ham.get_mei(occ_to_site[i], occ_to_site[i]);

          // (ii|ii) and (ii|jj)
          for (int j = 0; j < nap; j++)
            retval += 0.5 * rat * ham.get_tei(occ_to_site[i], occ_to_site[i], occ_to_site[j], occ_to_site[j]);

          // (ia|ai)
          for (int a = 0; a < nau; a++)
            retval += 0.5 * rat * ham.get_tei(occ_to_site[i], unocc_to_site[a], unocc_to_site[a], occ_to_site[i]);

        }

      }

      //formic::start_timer("same spin single and double excite");

      const bool slow_aaaa = false;
      //const bool slow_aaaa = true;

      // slow version
      if ( slow_aaaa || ( der_coeff != 0 && cont_inv != 0 ) ) {

        // off-diagonal one electron terms and half-diag two electron terms
        for (int i = 0; i < nap; i++) {

          // get site where we destroy an electron
          newly_empt[0] = 2*occ_to_site[i] + beta;

          for (int a = 0; a < nau; a++) {

            // get site where we create an electron
            newly_full[0] = 2*unocc_to_site[a] + beta;

            // get coefficient ratio
            const S rat = formic::fqmc::rgp_det_ratio_deriv(nap,
                                                            nas,
                                                            2,
                                                            newly_full,
                                                            newly_empt,
                                                            gamma,
                                                            inv_mat,
                                                            der_coeff,
                                                            cont_inv,
                                                            a_occ_to_site,
                                                            b_occ_to_site,
                                                            a_site_to_occ,
                                                            b_site_to_occ,
                                                            work);

            // ia
            retval += rat * ham.get_mei(occ_to_site[i], unocc_to_site[a]);

            // (ii|ia) and (ii|ja), although the code exchanges i and j
            for (int j = 0; j < nap; j++)
              retval += 0.5 * rat * ham.get_tei(occ_to_site[j], occ_to_site[j], occ_to_site[i], unocc_to_site[a]);

            // (ia|jj)
            for (int j = 0; j < nap; j++)
              retval += 0.5 * rat * ham.get_tei(occ_to_site[i], unocc_to_site[a], occ_to_site[j], occ_to_site[j]);
            retval -= 0.5 * rat * ham.get_tei(occ_to_site[i], unocc_to_site[a], occ_to_site[i], occ_to_site[i]); // remove (ia|ii)

            // (ia|aa) and (ia|ab), although the code exchanges a and b
            for (int b = 0; b < nau; b++)
              retval += 0.5 * rat * ham.get_tei(occ_to_site[i], unocc_to_site[b], unocc_to_site[b], unocc_to_site[a]);

            // (ia|ji), although the code exchanges i and j
            for (int j = 0; j < nap; j++)
              retval -= 0.5 * rat * ham.get_tei(occ_to_site[j], unocc_to_site[a], occ_to_site[i], occ_to_site[j]);
            retval += 0.5 * rat * ham.get_tei(occ_to_site[i], unocc_to_site[a], occ_to_site[i], occ_to_site[i]); // remove (ia|ii)

            // bbaa (jj|ia) and aabb (ia|jj)
            if (beta) {
              for (int j = 0; j < nap; j++)
                retval += rat * ham.get_tei(a_occ_to_site[j],a_occ_to_site[j],occ_to_site[i],unocc_to_site[a]);
            } else {
              for (int j = 0; j < nap; j++)
                retval += rat * ham.get_tei(b_occ_to_site[j],b_occ_to_site[j],occ_to_site[i],unocc_to_site[a]);
            }

          }

        }

        // off-diagonal two electron terms
        for (int i = 0; i < nap; i++) {
          for (int a = 0; a < nau; a++) {
            for (int j = 0; j < i; j++) {
              for (int b = 0; b < a; b++) {

                // get sites where we destroy an electron
                newly_empt[0] = 2*occ_to_site[i] + beta;
                newly_empt[1] = 2*occ_to_site[j] + beta;

                // get sites where we create an electron
                newly_full[0] = 2*unocc_to_site[a] + beta;
                newly_full[1] = 2*unocc_to_site[b] + beta;

                // get coefficient ratio
                const S rat = formic::fqmc::rgp_det_ratio_deriv(nap,
                                                                nas,
                                                                4,
                                                                newly_full,
                                                                newly_empt,
                                                                gamma,
                                                                inv_mat,
                                                                der_coeff,
                                                                cont_inv,
                                                                a_occ_to_site,
                                                                b_occ_to_site,
                                                                a_site_to_occ,
                                                                b_site_to_occ,
                                                                work);

                // (ia|jb)
                retval += 2.0 * 0.5 * rat * ham.get_tei(occ_to_site[i], unocc_to_site[a], occ_to_site[j], unocc_to_site[b]);
                retval -= 2.0 * 0.5 * rat * ham.get_tei(occ_to_site[i], unocc_to_site[b], occ_to_site[j], unocc_to_site[a]);

              }
            }
          }
        }

      // fast version
      } else {

        // get some arrays we will need
        const int no = nap;
        const int nv = nau;
//        static std::vector<S>    tia(no*nv, formic::zero(S()));
//        static std::vector<S>     ti(no, formic::zero(S()));
//        static std::vector<S>  viajb(no*nv*no*nv, formic::zero(S()));
//        static std::vector<S>   vijb(no*no*nv, formic::zero(S()));
//        static std::vector<S>    vij(no*no, formic::zero(S()));
//        static std::vector<S>    via(no*nv, formic::zero(S()));
//        static std::vector<S>     vi(no, formic::zero(S()));
//        static std::vector<S>   Vmat(nv*no, formic::zero(S()));
//        static std::vector<S>   Wmat(no*no, formic::zero(S()));
//        static std::vector<S>   Bmat(no*no, formic::zero(S()));
//        static std::vector<S>  VBmat(no*nv, formic::zero(S()));
//        static std::vector<S>  WBmat(no*no, formic::zero(S()));
//        static std::vector<S>    wia(no*nv, formic::zero(S()));
//        static std::vector<S>    wij(no*no, formic::zero(S()));
//        static std::vector<S>    waj(nv*no, formic::zero(S()));

        static S * big_array = new S[ no*nv
                                     +no
                                     +no*nv*no*nv
                                     +no*no*nv
                                     +no*no
                                     +no*nv
                                     +no
                                     +nv*no
                                     +no*no
                                     +no*no
                                     +no*nv
                                     +no*no
                                     +no*nv
                                     +no*no
                                     +nv*no ];
        S * next_pos = big_array;
        static S *    tia = next_pos;    next_pos += no*nv;
        static S *     ti = next_pos;    next_pos += no;
        static S *  viajb = next_pos;    next_pos += no*nv*no*nv;
        static S *   vijb = next_pos;    next_pos += no*no*nv;
        static S *    vij = next_pos;    next_pos += no*no;
        static S *    via = next_pos;    next_pos += no*nv;
        static S *     vi = next_pos;    next_pos += no;
        static S *   Vmat = next_pos;    next_pos += nv*no;
        static S *   Wmat = next_pos;    next_pos += no*no;
        static S *   Bmat = next_pos;    next_pos += no*no;
        static S *  VBmat = next_pos;    next_pos += no*nv;
        static S *  WBmat = next_pos;    next_pos += no*no;
        static S *    wia = next_pos;    next_pos += no*nv;
        static S *    wij = next_pos;    next_pos += no*no;
        static S *    waj = next_pos;    next_pos += nv*no;

        static S * full_tei = new S[ nas*nas*nas*nas ];

        static bool to_initialize = true;
        if (to_initialize) {

          std::fill(big_array, big_array + no*nv
                                         + no
                                         + no*nv*no*nv
                                         + no*no*nv
                                         + no*no
                                         + no*nv
                                         + no
                                         + nv*no
                                         + no*no
                                         + no*no
                                         + no*nv
                                         + no*no
                                         + no*nv
                                         + no*no
                                         + nv*no,          formic::zero(S()));

          for (int a = 0; a < nas; a++)
          for (int b = 0; b <= a ; b++)
          for (int c = 0; c < nas; c++)
          for (int d = 0; d <= c ; d++) {
            const S integ = ham.get_tei(a,b,c,d);
            full_tei[a+nas*(b+nas*(c+nas*d))] = integ;
            full_tei[a+nas*(b+nas*(d+nas*c))] = integ;
            full_tei[b+nas*(a+nas*(c+nas*d))] = integ;
            full_tei[b+nas*(a+nas*(d+nas*c))] = integ;
          }

          to_initialize = false;

        }

        // prepare inverse matrix
        if (beta) {
          for (int j = 0; j < nap; j++)
          for (int i = 0; i < nap; i++)
            Bmat[i+nap*j] = inv_mat[j*nap+i];
        } else {
          for (int j = 0; j < nap; j++)
          for (int i = 0; i < nap; i++)
            Bmat[i+nap*j] = inv_mat[i*nap+j];
        }

        // prepare new row matrix
        if (beta) {
          for (int a = 0; a < nau; a++)
          for (int i = 0; i < nap; i++)
            Vmat[i+nap*a] = gamma[ unocc_to_site[a] + nas * a_occ_to_site[i] ];
        } else {
          for (int a = 0; a < nau; a++)
          for (int i = 0; i < nap; i++)
            Vmat[i+nap*a] = gamma[ unocc_to_site[a] * nas + b_occ_to_site[i] ];
        }

        // prepare old row matrix
        if (beta) {
          for (int j = 0; j < nap; j++)
          for (int i = 0; i < nap; i++)
            Wmat[i+nap*j] = gamma[ occ_to_site[j] + nas * a_occ_to_site[i] ];
        } else {
          for (int j = 0; j < nap; j++)
          for (int i = 0; i < nap; i++)
            Wmat[i+nap*j] = gamma[ occ_to_site[j] * nas + b_occ_to_site[i] ];
        }

//        // prepare one electron integrals
//        for (int a = 0; a < nau; a++)
//        for (int i = 0; i < nap; i++) {
//
//          // start with zero
//          tia[i+nap*a] = formic::zero(S());
//
//          // ia one electron integral
//          tia[i+nap*a] += ham.get_mei(occ_to_site[i], unocc_to_site[a]);
//
//          for (int j = 0; j < nap; j++) {
//
//            const S integral = full_tei[occ_to_site[j]+nas*(occ_to_site[j]+nas*(occ_to_site[i]+nas*unocc_to_site[a]))];
//            //const S integral = ham.get_tei(occ_to_site[j], occ_to_site[j], occ_to_site[i], unocc_to_site[a]);
//
//            // (ii|ia) and (ii|ja), although the code exchanges i and j
//            tia[i+nap*a] += 0.5 * integral;
//
//            // (ia|jj)
//            tia[i+nap*a] += 0.5 * integral;
//
//          }
//
//          //tia[i+nap*a] -= 0.5 * ham.get_tei(occ_to_site[i], unocc_to_site[a], occ_to_site[i], occ_to_site[i]); // remove (ia|ii)
//
//          // (ia|aa) and (ia|ab), although the code exchanges a and b
//          for (int b = 0; b < nau; b++)
//            tia[i+nap*a] += 0.5 * full_tei[occ_to_site[i]+nas*(unocc_to_site[b]+nas*(unocc_to_site[b]+nas*unocc_to_site[a]))];
//            //tia[i+nap*a] += 0.5 * ham.get_tei(occ_to_site[i], unocc_to_site[b], unocc_to_site[b], unocc_to_site[a]);
//
//          // (ia|ji), although the code exchanges i and j
//          for (int j = 0; j < nap; j++)
//            tia[i+nap*a] -= 0.5 * full_tei[occ_to_site[j]+nas*(unocc_to_site[a]+nas*(occ_to_site[i]+nas*occ_to_site[j]))];
//            //tia[i+nap*a] -= 0.5 * ham.get_tei(occ_to_site[j], unocc_to_site[a], occ_to_site[i], occ_to_site[j]);
//
//          //tia[i+nap*a] += 0.5 * ham.get_tei(occ_to_site[i], unocc_to_site[a], occ_to_site[i], occ_to_site[i]); // remove (ia|ii)
//
//        }

        // prepare one electron integrals

        // ia one electron integral
        for (int a = 0; a < nau; a++)
        for (int i = 0; i < nap; i++)
          tia[i+nap*a] = ham.get_mei(occ_to_site[i], unocc_to_site[a]);

        for (int a = 0; a < nau; a++)
        for (int i = 0; i < nap; i++)
        for (int j = 0; j < nap; j++) {

          const S integral = full_tei[occ_to_site[j]+nas*(occ_to_site[j]+nas*(occ_to_site[i]+nas*unocc_to_site[a]))];
          //const S integral = ham.get_tei(occ_to_site[j], occ_to_site[j], occ_to_site[i], unocc_to_site[a]);

          // (ii|ia) and (ii|ja), although the code exchanges i and j
          tia[i+nap*a] += 0.5 * integral;

          // (ia|jj)
          tia[i+nap*a] += 0.5 * integral;

        }

        // (ia|aa) and (ia|ab), although the code exchanges a and b
        for (int a = 0; a < nau; a++)
        for (int b = 0; b < nau; b++)
        for (int i = 0; i < nap; i++)
          tia[i+nap*a] += 0.5 * full_tei[occ_to_site[i]+nas*(unocc_to_site[b]+nas*(unocc_to_site[b]+nas*unocc_to_site[a]))];
          //tia[i+nap*a] += 0.5 * full_tei[occ_to_site[i]+nas*(unocc_to_site[b]+nas*(unocc_to_site[a]+nas*unocc_to_site[b]))];

        // (ia|ji), although the code exchanges i and j
        for (int a = 0; a < nau; a++)
        for (int j = 0; j < nap; j++)
        for (int i = 0; i < nap; i++)
          tia[i+nap*a] -= 0.5 * full_tei[occ_to_site[i]+nas*(occ_to_site[j]+nas*(occ_to_site[j]+nas*unocc_to_site[a]))];
          //tia[i+nap*a] -= 0.5 * full_tei[occ_to_site[j]+nas*(unocc_to_site[a]+nas*(occ_to_site[i]+nas*occ_to_site[j]))];

        // bbaa (jj|ia) and aabb (ia|jj)
        if (beta) {
          for (int a = 0; a < nau; a++)
          for (int i = 0; i < nap; i++)
          for (int j = 0; j < nap; j++)
            tia[i+nap*a] += full_tei[a_occ_to_site[j]+nas*(a_occ_to_site[j]+nas*(occ_to_site[i]+nas*unocc_to_site[a]))];
        } else {
          for (int a = 0; a < nau; a++)
          for (int i = 0; i < nap; i++)
          for (int j = 0; j < nap; j++)
            tia[i+nap*a] += full_tei[b_occ_to_site[j]+nas*(b_occ_to_site[j]+nas*(occ_to_site[i]+nas*unocc_to_site[a]))];
        }

        // prepare two electron integrals
        std::fill(viajb, viajb + no*nv*no*nv, formic::zero(S()));
        for (int b = 0; b < nau; b++)
        for (int j = 0; j < nap; j++)
        for (int a = 0; a < b  ; a++)
        for (int i = 0; i < j  ; i++) {
          const S integ1 = full_tei[occ_to_site[i]+nas*(unocc_to_site[a]+nas*(occ_to_site[j]+nas*unocc_to_site[b]))];
          const S integ2 = full_tei[occ_to_site[j]+nas*(unocc_to_site[a]+nas*(occ_to_site[i]+nas*unocc_to_site[b]))];
          //const S integ1 = ham.get_tei(occ_to_site[i], unocc_to_site[a], occ_to_site[j], unocc_to_site[b]);
          //const S integ2 = ham.get_tei(occ_to_site[j], unocc_to_site[a], occ_to_site[i], unocc_to_site[b]);
          viajb[i+nap*(a+nau*(j+nap*b))] = integ1;
          viajb[j+nap*(b+nau*(i+nap*a))] = integ1;
          viajb[j+nap*(a+nau*(i+nap*b))] = integ2;
          viajb[i+nap*(b+nau*(j+nap*a))] = integ2;
        }

        //// prepare intermediates
        //prep_aaaa(
        //          no,
        //          nv,
        //          &tia[0],
        //          &ti[0],
        //          &viajb[0],
        //          &vijb[0],
        //          &vij[0],
        //          &via[0],
        //          &vi[0],
        //          &Vmat[0],
        //          &Wmat[0],
        //          &Bmat[0],
        //          &VBmat[0],
        //          &WBmat[0]);

        // compute contribution from single and double excitations
        get_aaaa(
                 no,
                 nv,
                 &tia[0],
                 //&ti[0],
                 &viajb[0],
                 //&vijb[0],
                 //&vij[0],
                 //&via[0],
                 //&vi[0],
                 &Bmat[0],
                 &Vmat[0],
                 &VBmat[0],
                 //&WBmat[0],
                 //&wia[0],
                 //&wij[0],
                 &waj[0],
                 &retval);

      }

      //formic::stop_timer("same spin single and double excite");

      // return the contribution
      return retval;

    }

    // function to compute energy from opposite-spin Hamiltonian operators
    static S opp_spin(const int beta,
                      const int nap,
                      const int nas,
                      const formic::Hamiltonian<S> & ham,
                      const S * const gamma,
                      const S * const inv_mat,
                      const S * const der_coeff,
                      const S * const cont_inv,
                      const int * const a_occ_to_site,
                      const int * const b_occ_to_site,
                      const int * const a_site_to_occ,
                      const int * const b_site_to_occ,
                      const int * const a_unocc_to_site,
                      const int * const b_unocc_to_site,
                      S * const work)
    {

      const int nau = nas - nap;

      // get work vectors
      S * const wv1 = work + 0*nap;
      S * const wv2 = work + 1*nap;
      S * const wv3 = work + 2*nap;

      // initialize contribution
      S retval = formic::zero(S());

      // get integer arrays describe the excitations
      int newly_full[] = {0,0};
      int newly_empt[] = {0,0};

      // diagonal two electron terms
      if ( der_coeff == 0 || cont_inv == 0 ) {

        S rat;
        if ( beta == 0 )
          rat = formic::fqmc::rgp_det_ratio_deriv(nap,
                                                  nas,
                                                  0,
                                                  newly_full,
                                                  newly_empt,
                                                  gamma,
                                                  inv_mat,
                                                  der_coeff,
                                                  cont_inv,
                                                  a_occ_to_site,
                                                  b_occ_to_site,
                                                  a_site_to_occ,
                                                  b_site_to_occ,
                                                  work);
        else
          rat = formic::fqmc::rgp_det_ratio_deriv(nap,
                                                  nas,
                                                  0,
                                                  newly_full,
                                                  newly_empt,
                                                  gamma,
                                                  inv_mat,
                                                  der_coeff,
                                                  cont_inv,
                                                  b_occ_to_site,
                                                  a_occ_to_site,
                                                  b_site_to_occ,
                                                  a_site_to_occ,
                                                  work);

        // (ii|jj)
        for (int i = 0; i < nap; i++)
          for (int j = 0; j < nap; j++)
            retval += 0.5 * rat * ham.get_tei(a_occ_to_site[i], a_occ_to_site[i], b_occ_to_site[j], b_occ_to_site[j]);

      }

      //formic::start_timer("opp spin single and double excite");

//      // diag alpha, off-diag beta
//      if ( der_coeff != 0 && cont_inv != 0 ) {
//        for (int i = 0; i < nap; i++) {
//          for (int a = 0; a < nau; a++) {
//
//            // get site where we destroy an electron
//            newly_empt[0] = 2*b_occ_to_site[i] + 1-beta;
//
//            // get site where we create an electron
//            newly_full[0] = 2*b_unocc_to_site[a] + 1-beta;
//
//            // get coefficient ratio
//            S rat;
//            if ( beta == 0 )
//              rat = formic::fqmc::rgp_det_ratio_deriv(nap,
//                                                      nas,
//                                                      2,
//                                                      newly_full,
//                                                      newly_empt,
//                                                      gamma,
//                                                      inv_mat,
//                                                      der_coeff,
//                                                      cont_inv,
//                                                      a_occ_to_site,
//                                                      b_occ_to_site,
//                                                      a_site_to_occ,
//                                                      b_site_to_occ,
//                                                      work);
//            else
//              rat = formic::fqmc::rgp_det_ratio_deriv(nap,
//                                                      nas,
//                                                      2,
//                                                      newly_full,
//                                                      newly_empt,
//                                                      gamma,
//                                                      inv_mat,
//                                                      der_coeff,
//                                                      cont_inv,
//                                                      b_occ_to_site,
//                                                      a_occ_to_site,
//                                                      b_site_to_occ,
//                                                      a_site_to_occ,
//                                                      work);
//
//            // (jj|ia) and (ia|jj)
//            for (int j = 0; j < nap; j++)
//              retval += 2.0 * 0.5 * rat * ham.get_tei(a_occ_to_site[j], a_occ_to_site[j], b_occ_to_site[i], b_unocc_to_site[a]);
//
//          }
//        }
//      }

      const bool slow_aabb = false;
      //const bool slow_aabb = true;

      // off-diag alpha, off-diag beta
      if ( beta == 0 && ( slow_aabb || ( der_coeff != 0 && cont_inv != 0 ) ) ) {

        for (int i = 0; i < nap; i++) {
          for (int a = 0; a < nau; a++) {
            for (int j = 0; j < nap; j++) {
              for (int b = 0; b < nau; b++) {

                // get sites where we destroy an electron
                newly_empt[0] = 2*a_occ_to_site[i] + beta;
                newly_empt[1] = 2*b_occ_to_site[j] + 1-beta;

                // get sites where we create an electron
                newly_full[0] = 2*a_unocc_to_site[a] + beta;
                newly_full[1] = 2*b_unocc_to_site[b] + 1-beta;

                // get coefficient ratio
                const S rat = formic::fqmc::rgp_det_ratio_deriv(nap,
                                                                nas,
                                                                4,
                                                                newly_full,
                                                                newly_empt,
                                                                gamma,
                                                                inv_mat,
                                                                der_coeff,
                                                                cont_inv,
                                                                a_occ_to_site,
                                                                b_occ_to_site,
                                                                a_site_to_occ,
                                                                b_site_to_occ,
                                                                work);

                // (ia|jb)
                retval += 2.0 * 0.5 * rat * ham.get_tei(a_occ_to_site[i], a_unocc_to_site[a], b_occ_to_site[j], b_unocc_to_site[b]);

              }
            }
          }
        }

      } else if ( beta == 0 ) {

        const int no = nap;
        const int nv = nau;

        static S * big_array = new S[  no*nv*no*nv
                                     + no*nv
                                     + no*nv
                                     + nv*nv
                                     + no*no
                                     + no*no
                                     + no*nv
                                     + no*nv
                                     + nv*nv
                                     + nv*nv
                                    ];

        S * next_pos = big_array;
        static S *  Tiajb = next_pos; next_pos += no*nv*no*nv;
        static S * newrow = next_pos; next_pos += no*nv;
        static S * newcol = next_pos; next_pos += no*nv;
        static S *   Xmat = next_pos; next_pos += nv*nv;
        static S *   Bmat = next_pos; next_pos += no*no;
        static S *   Btrn = next_pos; next_pos += no*no;
        static S *   VBia = next_pos; next_pos += no*nv;
        static S *   BUkb = next_pos; next_pos += no*nv;
        static S *  VBUab = next_pos; next_pos += nv*nv;
        static S *   TBab = next_pos; next_pos += nv*nv;

        static S * full_tei = new S[ nas*nas*nas*nas ];

        static bool to_initialize = true;
        if (to_initialize) {

          std::fill(big_array, big_array + no*nv*no*nv
                                         + no*nv
                                         + no*nv
                                         + nv*nv
                                         + no*no
                                         + no*no
                                         + no*nv
                                         + no*nv
                                         + nv*nv
                                         + nv*nv, formic::zero(S()));

          for (int a = 0; a < nas; a++)
          for (int b = 0; b <= a ; b++)
          for (int c = 0; c < nas; c++)
          for (int d = 0; d <= c ; d++) {
            const S integ = ham.get_tei(a,b,c,d);
            full_tei[a+nas*(b+nas*(c+nas*d))] = integ;
            full_tei[a+nas*(b+nas*(d+nas*c))] = integ;
            full_tei[b+nas*(a+nas*(c+nas*d))] = integ;
            full_tei[b+nas*(a+nas*(d+nas*c))] = integ;
          }

          to_initialize = false;

        }

        // prepare inverse matrix
        for (int j = 0; j < nap; j++)
        for (int i = 0; i < nap; i++)
          Bmat[i+nap*j] = inv_mat[i*nap+j];

        // prepare new row matrix
        for (int a = 0; a < nau; a++)
        for (int i = 0; i < nap; i++)
          newrow[i+nap*a] = gamma[ a_unocc_to_site[a] * nas + b_occ_to_site[i] ];

        // prepare new column matrix
        for (int a = 0; a < nau; a++)
        for (int i = 0; i < nap; i++)
          newcol[i+nap*a] = gamma[ a_occ_to_site[i] * nas + b_unocc_to_site[a] ];

        // prepare unoccupied pairing matrix
        for (int b = 0; b < nau; b++)
        for (int a = 0; a < nau; a++)
          Xmat[a+nau*b] = gamma[ a_unocc_to_site[a] * nas + b_unocc_to_site[b] ];

        // prepare two electron integrals
        for (int b = 0; b < nau; b++)
        for (int j = 0; j < nap; j++)
        for (int a = 0; a < nau; a++)
        for (int i = 0; i < nap; i++) {
          const S integ = full_tei[a_occ_to_site[i]+nas*(a_unocc_to_site[a]+nas*(b_occ_to_site[j]+nas*b_unocc_to_site[b]))];
          Tiajb[i+nap*(a+nau*(j+nap*b))] = integ;
        }

        // get energy contribution
        get_aabb(
                 no,
                 nv,
                &Tiajb[0],
               &newrow[0],
               &newcol[0],
                 &Xmat[0],
                 &Bmat[0],
                 &Btrn[0],
                 &VBia[0],
                 &BUkb[0],
                &VBUab[0],
                 &TBab[0],
                 &retval);

      }

      //formic::stop_timer("opp spin single and double excite");

      // return the contribution
      return retval;

    }

  };

  // initialize local energy
  S retval = formic::zero(S());

  // compute same-spin terms
  retval += local_funcs::same_spin(0, nap, nas, ham, gamma, inv_mat, der_coeff, cont_inv, a_occ_to_site, a_unocc_to_site,
                                   a_occ_to_site, b_occ_to_site, a_site_to_occ, b_site_to_occ, work);
  retval += local_funcs::same_spin(1, nap, nas, ham, gamma, inv_mat, der_coeff, cont_inv, b_occ_to_site, b_unocc_to_site,
                                   a_occ_to_site, b_occ_to_site, a_site_to_occ, b_site_to_occ, work);

  // compute opposite-spin terms
  retval += local_funcs::opp_spin(0, nap, nas, ham, gamma, inv_mat, der_coeff, cont_inv,
                                  a_occ_to_site, b_occ_to_site,
                                  a_site_to_occ, b_site_to_occ,
                                  a_unocc_to_site, b_unocc_to_site,
                                  work);
  retval += local_funcs::opp_spin(1, nap, nas, ham, gamma, inv_mat, der_coeff, cont_inv,
                                  b_occ_to_site, a_occ_to_site,
                                  b_site_to_occ, a_site_to_occ,
                                  b_unocc_to_site, a_unocc_to_site,
                                  work);

  // add nuclear repulsion term
  if ( der_coeff == 0 || cont_inv == 0 )
    retval += ham.e_cor_nuc() * formic::fqmc::rgp_det_ratio_deriv(nap,
                                                                  nas,
                                                                  0,
                                                                  (int *)0,
                                                                  (int *)0,
                                                                  gamma,
                                                                  inv_mat,
                                                                  der_coeff,
                                                                  cont_inv,
                                                                  a_occ_to_site,
                                                                  b_occ_to_site,
                                                                  a_site_to_occ,
                                                                  b_site_to_occ,
                                                                  work);

  // return the local energy (or Hamiltonian-based derivative contraction)
  return retval;

}

template double formic::fqmc::rgp_simple_local_energy(const int nap,
                                                      const int nas,
                                                      const formic::Hamiltonian<double> & ham,
                                                      const double * const gamma,
                                                      const double * const inv_mat,
                                                      const double * const der_coeff,
                                                      const double * const cont_inv,
                                                      const int * const a_occ_to_site,
                                                      const int * const b_occ_to_site,
                                                      const int * const a_site_to_occ,
                                                      const int * const b_site_to_occ,
                                                      const int * const a_unocc_to_site,
                                                      const int * const b_unocc_to_site,
                                                      const int * const a_site_to_unocc,
                                                      const int * const b_site_to_unocc,
                                                      double * const work);
template std::complex<double>  formic::fqmc::rgp_simple_local_energy(const int nap,
                                                                     const int nas,
                                                                     const formic::Hamiltonian<std::complex<double> > & ham,
                                                                     const std::complex<double>  * const gamma,
                                                                     const std::complex<double>  * const inv_mat,
                                                                     const std::complex<double> * const der_coeff,
                                                                     const std::complex<double> * const cont_inv,
                                                                     const int * const a_occ_to_site,
                                                                     const int * const b_occ_to_site,
                                                                     const int * const a_site_to_occ,
                                                                     const int * const b_site_to_occ,
                                                                     const int * const a_unocc_to_site,
                                                                     const int * const b_unocc_to_site,
                                                                     const int * const a_site_to_unocc,
                                                                     const int * const b_site_to_unocc,
                                                                     std::complex<double>  * const work);
